﻿'
' Copyright (c) 2009 Aviation Data Systems, Inc
' http://fboweb.com
'
' This software is licensed under the GNU Public License v2
'
'
'
'


Public Class point

    ''' <summary>
    ''' The point's latitude.
    ''' </summary>
    Private _lat As Double
    ''' <summary>
    ''' The point's longitude.
    ''' </summary>
    Private _lon As Double

    ''' <summary>
    ''' The point's altitude (in feet).
    ''' </summary>
    Private _alt As Int32

    Public Sub New(ByVal lat As Double, ByVal lon As Double)
        _lat = lat
        _lon = lon
    End Sub

    ''' <summary>
    ''' Gets or sets the latitude of the point.
    ''' </summary>
    Public Property Lat() As Double
        Get
            Return _lat
        End Get
        Set(ByVal value As Double)
            _lat = value
        End Set
    End Property

    ''' <summary>
    ''' Gets or sets the longitude of the point.
    ''' </summary>
    Public Property Lon() As Double
        Get
            Return _lon
        End Get
        Set(ByVal value As Double)
            _lon = value
        End Set
    End Property

    ''' <summary>
    ''' Gets or sets the point's altitude (in feet).
    ''' </summary>
    Public Property Altitude() As Int32
        Get
            Return _alt
        End Get
        Set(ByVal value As Int32)
            _alt = value
        End Set
    End Property

    ''' <summary>
    ''' Returns the true course to the specified point.
    ''' </summary>
    Public Function CourseTo(ByVal point As point) As Short
        Return CourseTo(point.Lat, point.Lon)
    End Function

    ''' <summary>
    ''' Returns the true course to the specified point.
    ''' </summary>
    Public Function CourseTo(ByVal lat As Double, ByVal lon As Double) As Short
        If lat = Me.Lat AndAlso lon = Me.Lon Then
            Return 0
        Else
            Try
                Dim lat1, long1, lat2, long2 As Double

                lat1 = Me.Lat / 60
                long1 = Me.Lon / 60
                lat2 = lat / 60
                long2 = lon / 60

                Dim dlong, dlat, a, d, heading As Double
                dlong = long2 - long1
                dlat = lat2 - lat1

                a = Math.Sin(dlat / 2) ^ 2 + Math.Cos(lat1) * Math.Cos(lat2) * Math.Sin(dlong / 2) ^ 2
                d = 2 * Math.Atan(Math.Sqrt(a) / Math.Sqrt(1 - a))

                heading = Math.Acos((Math.Sin(lat2) - Math.Sin(lat1) * Math.Cos(d)) / (Math.Sin(d) * Math.Cos(lat1)))

                If heading = heading Then ' a little trick, when heading=-1.#IND or whatever, this is false.  This traps coordinates that are too close to plot.
                    If Math.Sin(long2 - long1) < 0 Then
                        heading = 2 * 3.1415926535897931 - heading
                    End If

                    a = radians_to_degrees(heading)
                    If a < 0 Then a = 180 + a
                    If a > 360 Then a = a - 180

                    CourseTo = a
                Else
                    ' too close to compute
                    CourseTo = 0
                End If
            Catch ex As Exception
                ' too close to compute
                CourseTo = 0
            End Try
        End If
    End Function

    Function radians_to_degrees(ByVal a)
        radians_to_degrees = (180 / 3.1415926535897931) * a
    End Function

    ''' <summary>
    ''' Returns the distance (in nautical miles) between two points.
    ''' </summary>
    Public Function Distance(ByVal latpoint1 As Double, ByVal longpoint1 As Double, ByVal latpoint2 As Double, ByVal longpoint2 As Double) As Double
        Dim meters_per_mile As Double = 1609.344
        Dim nautical_miles_per_mile As Double = (5280 / 6076)
        Dim la1, lo1, la2, lo2 As Double

        la1 = degrees_to_radians(latpoint1)
        lo1 = degrees_to_radians(longpoint1)
        la2 = degrees_to_radians(latpoint2)
        lo2 = degrees_to_radians(longpoint2)

        Distance = (great_circle_distance(la1, lo1, la2, lo2) / meters_per_mile) * nautical_miles_per_mile
    End Function

    Private Function great_circle_distance(ByVal lat1, ByVal long1, ByVal lat2, ByVal long2) As Double
        '//	' Approx radius of Earth in meters.  True radius varies from 6357km (polar) to 6378km (equatorial)
        Dim earth_radius As Double = 6367000.0
        Dim dlong, dlat, a, d As Double

        dlong = long2 - long1
        dlat = lat2 - lat1

        a = Math.Sin(dlat / 2) ^ 2 + Math.Cos(lat1) * Math.Cos(lat2) * Math.Sin(dlong / 2) ^ 2
        d = 2 * Math.Atan(Math.Sqrt(a) / Math.Sqrt(1 - a))
        '//	d=2*asin(sqr((Math.sin((dlat)/2))^2 + Math.cos(lat1)*Math.cos(lat2)*(Math.sin((dlong)/2))^2));

        great_circle_distance = earth_radius * d
    End Function

    ''' <summary>
    ''' Returns the distance (in nautical miles) from the specified point.
    ''' </summary>
    Public Function DistanceFrom(ByVal p As point) As Double
        DistanceFrom = Distance(Me.Lat, Me.Lon, p.Lat, p.Lon)
    End Function

    ''' <summary>
    ''' Returns the distance (in nautical miles) from the specified point.
    ''' </summary>
    Public Function DistanceFrom(ByVal lat As Double, ByVal lon As Double)
        DistanceFrom = Distance(Me.Lat, Me.Lon, lat, lon)
    End Function

    Function degrees_to_radians(ByVal a)
        degrees_to_radians = a * (3.1415926535897931 / 180.0)
    End Function
End Class
